library(lmerTest)
library(psych)
library(ggplot2)
library(simpleboot)
library(boot)
library(xtable)
library(brms)
source("BrmsTables.R")
test <- read.csv("CovidExp.csv")
completed <- test$participant._index_in_pages >= 20

test <- subset(test, completed)
nrow(test)
test$ageHigh <- test$player.age >= median(test$player.age,na.rm = T) 
test$educationHigh <- test$player.education == "Laurea o superiore" 
test$female <- ifelse(test$player.sex == "Femmina",1,0)
test$player.covidThreat <- abs(test$player.covidThreat - 9)
test$player.orderDonation <- as.factor(test$player.orderDonation)
test$donateDum <- test$player.donate != 1
test$pol <- test$player.treatment == 2 | test$player.treatment == 4
test$sci <- test$player.treatment == 3 | test$player.treatment == 4

#reducing many questions to few factors
fa(test[,38:44])$loadings
test$questions <- fa(test[,38:44])$scores
#reducing the three information requests to one factor
bla1 <- ifelse(test$player.moreInfoText == "Si",1,0)
bla2 <- ifelse(test$player.moreInfoInterventions == "Si",1,0)
bla3 <- ifelse(test$player.moreInfoStatements == "Si",1,0)
fa(cbind(bla1,bla2,bla3))$loadings
test$infoReq <- fa(cbind(bla1,bla2,bla3))$scores

#demographic subgroups for the analysis
women <-subset(test, test$player.sex == "Femmina")
men <- subset(test, test$player.sex == "Maschio")
old <- subset(test, test$ageHigh)
young <- subset(test, !test$ageHigh)

#demographics of the sample
library(latticeExtra)
demographicFilter <-  test$player.age >= 18
demo <- subset(test, demographicFilter)
table(demo$female)
library(plyr)
demo$player.employment <- revalue(demo$player.employment, c("Disoccupata/o"="Unemp", "Occupata/o autonoma/o"="Indep","Occupata/o dipendente privato" = "Priv",
                           "Occupata/o dipendente pubblico" = "Pub", "Pensionata/o" = "Pens", "Studente(ssa)" = "Stud"))
#figures S1-S3
histogram(~player.employment|player.treatment, data = demo,
          panel = function(x,...) {
            panel.histogram(x,..., col = "white")
          })
histogram(~demo$player.age,
          xlab = "age",
          panel = function(x,...) {
            panel.histogram(x,..., col = "white")
            panel.lines( c( median(x, na.rm = T)-0.5, median(x, na.rm = T)-0.5),c(-10, 100), col = "red", lwd = 3,cex = 0.1, lty = 2)
          })  #600 X 350
histogram(~I(demo$player.covidThreat + 1),
     xlab = "agreement", 
     main = "I perceive COVID-19 as a threat to my health",
     breaks = 0:10,
     panel = function(x,...) {
       panel.histogram(x,..., col = "white")
       panel.lines( c( median(x, na.rm = T) - 0.5, median(x, na.rm = T) - 0.5),c(-10, 100), col = "red", lwd = 3,cex = 0.1, lty = 2)
     }) #600 X 350
table(demo$player.covidFamiliy)
histogram(~demo$player.pol,
          xlab = "left <--> right", 
          main = "Where do you see yourself on a left-right political spectrum",
          breaks = 0:9,
          panel = function(x,...) {
            panel.histogram(x,..., col = "white")
            panel.lines( c( median(x, na.rm = T)-0.5, median(x, na.rm = T)-0.5),c(-10, 100), col = "red", lwd = 3,cex = 0.1, lty = 2)
          })
table(demo$player.voted)

#demographic distribution across treatments
demos1 <- aggregate(list(demo$player.age,demo$player.covidThreat,demo$player.pol,
                         demo$female,demo$player.covidFamiliy == "Si",demo$player.voted == "Si",demo$educationHigh), 
                    by = list(demo$pol, demo$sci), FUN = mean, na.rm = T)
demos2 <- aggregate(list(demo$player.age,demo$player.covidThreat,demo$player.pol), by = list(demo$pol, demo$sci), FUN = sd, na.rm = T)
demos <- cbind(demos1[,1:3], demos2[,3],demos1[,4],demos2[,4],demos1[,5],demos2[,5], 
               demos1[,6], demos1[,7], demos1[,8], demos1[,9])
names(demos) <- c("Politicians","Scientists","age mean","sd","threat mean","sd","pol mean","sd",
                  "ratio female","family","voted","edu. high")
demos
#table S2
print(xtable(demos,floating=FALSE), include.rownames=F)
table(test$player.treatment)
table(test$donateDum)

#table S1
table(test$player.province)

##############################
#BAYES########################
##############################
library(brms)
seedBrm <- 666
iterations <- 10000
delta <- 0.99
modelDonate <- donateDum ~ pol * sci + (1|player.orderDonation)

donateBayes <- brm(modelDonate,
                      data = test,
                      family = bernoulli(link = "logit"), cores = 2, chains = 2,
                      control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
simple.tab(donateBayes)
donateBayesWomen <- brm(modelDonate,
                   data = women,
                   family = bernoulli(link = "logit"), cores = 2, chains = 2,
                   control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
donateBayesMen <- brm(modelDonate,
                   data = men,
                   family = bernoulli(link = "logit"), cores = 2, chains = 2,
                   control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
donateBayesYoung <- brm(modelDonate,
                   data = young,
                   family = bernoulli(link = "logit"), cores = 2, chains = 2,
                   control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
donateBayesOld <- brm(modelDonate,
                   data = old,
                   family = bernoulli(link = "logit"), cores = 2, chains = 2,
                   control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
donateComparison <- multi.tab(list(donateBayes, donateBayesWomen,donateBayesMen,donateBayesYoung,donateBayesOld), round.est=2)
names(donateComparison)[3:7] <- c("Overall","Women","Men","Young","Old")
#table 1
donateComparison


modelQuestions <- questions~pol * sci
questionsBayes <- brm(modelQuestions,
                   data = test,
                   cores = 2, chains = 2,
                   control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
simple.tab(questionsBayes)
questionsBayesWomen <- brm(modelQuestions,
                      data = women,
                      cores = 2, chains = 2,
                      control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
questionsBayesMen <- brm(modelQuestions,
                      data = men,
                      cores = 2, chains = 2,
                      control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
questionsBayesYoung <- brm(modelQuestions,
                      data = young,
                      cores = 2, chains = 2,
                      control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
questionsBayesOld <- brm(modelQuestions,
                      data = old,
                      cores = 2, chains = 2,
                      control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
questionsComparison <- multi.tab(list(questionsBayes, questionsBayesWomen,questionsBayesMen,questionsBayesYoung,questionsBayesOld), round.est=2)
names(questionsComparison)[3:7] <- c("Overall","Women","Men","Young","Old")
#table 3
questionsComparison
  
modelInfo <- infoReq ~ pol * sci
infoBayes <- brm(modelInfo,
                   data = test, 
                   cores = 2, chains = 2,
                   control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
simple.tab(infoBayes)
infoBayesWomen <- brm(modelInfo,
                 data = women, 
                 cores = 2, chains = 2,
                 control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
infoBayesMen <- brm(modelInfo,
                 data = men, 
                 cores = 2, chains = 2,
                 control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
infoBayesYoung <- brm(modelInfo,
                 data = young, 
                 cores = 2, chains = 2,
                 control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
infoBayesOld <- brm(modelInfo,
                 data = old, 
                 cores = 2, chains = 2,
                 control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
infoComparison <- multi.tab(list(infoBayes, infoBayesWomen,infoBayesMen,infoBayesYoung,infoBayesOld), round.est=2)
names(infoComparison)[3:7] <- c("Overall","Women","Men","Young","Old")
#table 2
infoComparison


#interaction plot 
library(gridExtra)
ci <- function(x, up=T){
  zz <- one.boot(x, mean, R=1000)
  i <- ifelse(up==T, 3, 2)
  boot.ci(zz, type="norm", conf=0.95)$normal[i]
}

dat <- test
dat$donation <- dat$donateDum
dat$politicians <- ifelse(dat$player.treatment %in% c(2,4), "Yes", "No")
dat$scientists <- ifelse(dat$player.treatment %in% c(3,4), "Yes", "No")
dat$opinions <- dat$questions

means <- aggregate(donation ~ politicians + scientists, dat, mean)
ci95u <- aggregate(donation ~ politicians + scientists, dat, ci)
ci95l <- aggregate(donation ~ politicians + scientists, dat, ci, up=F)
toplot1 <- data.frame(means[,1:2], mean=means[,3], ci95u=ci95u[,3], ci95l=ci95l[,3])
toplot1$variable <- "Donation proportion"

means <- aggregate(player.wantsInfo ~ politicians + scientists, dat, mean)
ci95u <- aggregate(player.wantsInfo ~ politicians + scientists, dat, ci)
ci95l <- aggregate(player.wantsInfo ~ politicians + scientists, dat, ci, up=F)
toplot2 <- data.frame(means[,1:2], mean=means[,3], ci95u=ci95u[,3], ci95l=ci95l[,3])
toplot2$variable <- "Info request proportion"

means <- aggregate(opinions ~ politicians + scientists, dat, mean)
ci95u <- aggregate(opinions ~ politicians + scientists, dat, ci)
ci95l <- aggregate(opinions ~ politicians + scientists, dat, ci, up=F)
toplot3 <- data.frame(means[,1:2], mean=means[,3], ci95u=ci95u[,3], ci95l=ci95l[,3])
toplot3$variable <- "Policy support"

toplot <- rbind(toplot1, toplot2, toplot3)
toplot$variable <- factor(toplot$variable, levels=c("Donation proportion","Info request proportion","Policy support"))

ggplot(toplot, aes(x=politicians, y=mean, color=scientists, group=scientists)) + 
  geom_point(size=4) + geom_line(size=1.5) + 
  geom_errorbar(aes(ymin=ci95l, ymax=ci95u), width=0.2) + 
  labs(color="Source of legitimacy: scientists", y=NULL, x="Source of legitimacy: politicians",
       caption="Whiskers show 95% CI") + facet_wrap(~ variable, scales="free") + 
  theme(legend.position="top", panel.grid=element_blank()) + scale_color_viridis_d() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.key=element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black")) +
   
#figure 2
ggsave("Interactions.png", width=25, height=12, units="cm")


##########Appendix####
#influence of the level of education
#first Bayesian

eduLow <- subset(test, !test$educationHigh)
eduHigh <- subset(test, test$educationHigh)

modelDonate <- donateDum ~ pol * sci + (1|player.orderDonation)
donateBayesEduLow <- brm(modelDonate,
                        data = eduLow,
                        family = bernoulli(link = "logit"), cores = 2, chains = 2,
                        control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
donateBayesEduHigh <- brm(modelDonate,
                      data = eduHigh,
                      family = bernoulli(link = "logit"), cores = 2, chains = 2,
                      control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)

donateComparisonEdu <- multi.tab(list(donateBayesEduLow,donateBayesEduHigh), round.est=2)
names(donateComparisonEdu)[3:4] <- c("lower educated","higher educated")
#table S4
donateComparisonEdu

modelInfo <- infoReq ~ pol * sci
infoBayesEduLow <- brm(modelInfo,
                         data = eduLow, cores = 2, chains = 2,
                         control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
infoBayesEduHigh <- brm(modelInfo,
                          data = eduHigh, cores = 2, chains = 2,
                          control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)

infoComparisonEdu <- multi.tab(list(infoBayesEduLow,infoBayesEduHigh), round.est=2)
names(infoComparisonEdu)[3:4] <- c("lower educated","higher educated")
infoComparisonEdu

modelQue <- questions ~ pol * sci
queBayesEduLow <- brm(modelQue,
                       data = eduLow, cores = 2, chains = 2,
                       control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
queBayesEduHigh <- brm(modelQue,
                        data = eduHigh, cores = 2, chains = 2,
                        control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)

queComparisonEdu <- multi.tab(list(queBayesEduLow,queBayesEduHigh), round.est=2)
names(queComparisonEdu)[3:4] <- c("lower educated","higher educated")
queComparisonEdu

totalEdu <- multi.tab(list(donateBayesEduLow,donateBayesEduHigh,infoBayesEduLow,infoBayesEduHigh,queBayesEduLow,queBayesEduHigh))
names(totalEdu)[3:8] <- c("lower educated","higher educated","lower educated","higher educated","lower educated","higher educated")
library(xtable)
print(xtable(totalEdu),floating=FALSE, include.rownames=FALSE)

#now the same but frequentist
donateFrqEduLow <- glmer(modelDonate,
                         data = eduLow,
                         family = "binomial")
donateFrqEduHigh <- glmer(modelDonate,
                          data = eduHigh,
                          family = "binomial")
infoFrqEduLow <- lm(modelInfo,
                       data = eduLow)
infoFrqEduHigh <- lm(modelInfo,
                        data = eduHigh)
queFrqEduLow <- lm(modelQue,
                      data = eduLow)
queFrwEduHigh <- lm(modelQue,
                       data = eduHigh)

stargazer(donateFrqEduLow,donateFrqEduHigh,infoFrqEduLow,infoFrqEduHigh,queFrqEduLow,queFrwEduHigh,
          star.cutoffs = c(0.05,0.01,0.001), model.numbers = FALSE,
          column.labels =  c("low edu","high edu","low edu","high edu","low edu","high edu"),
          omit.stat = c("adj.rsq", "f","ser","aic","bic","ll"),type="latex")
r2beta(donateFrqEduLow)
r2beta(donateFrqEduHigh)
r2beta(infoFrqEduLow)
r2beta(infoFrqEduHigh)
r2beta(queFrqEduLow)
r2beta(queFrwEduHigh)


#influence of the political stance
#first Bayesian
test$right <- test$player.pol >= median(test$player.pol, na.rm = T )
right <- subset(test, test$right)
left <- subset(test, !test$right)

modelDonate <- donateDum ~ pol * sci + (1|player.orderDonation)
donateBayesLeft <- brm(modelDonate,
                         data = left,
                         family = bernoulli(link = "logit"), cores = 2, chains = 2,
                         control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
donateBayesRight <- brm(modelDonate,
                          data = right,
                          family = bernoulli(link = "logit"), cores = 2, chains = 2,
                          control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)

donateComparisonPol <- multi.tab(list(donateBayesLeft,donateBayesRight), round.est=2)
names(donateComparisonPol)[3:4] <- c("left","right")
donateComparisonPol

modelInfo <- infoReq ~ pol * sci
infoBayesLeft <- brm(modelInfo,
                       data = left, cores = 2, chains = 2,
                       control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
infoBayesRight <- brm(modelInfo,
                        data = right, cores = 2, chains = 2,
                        control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)

infoComparisonPol <- multi.tab(list(infoBayesLeft,infoBayesRight), round.est=2)
names(infoComparisonPol)[3:4] <- c("left","right")
infoComparisonPol

modelQue <- questions ~ pol * sci
queBayesLeft <- brm(modelQue,
                      data = left, cores = 2, chains = 2,
                      control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)
queBayesRight <- brm(modelQue,
                       data = right, cores = 2, chains = 2,
                       control = list(adapt_delta = delta), iter = iterations, seed = seedBrm)

queComparisonPol <- multi.tab(list(queBayesLeft,queBayesRight), round.est=2)
names(queComparisonPol)[3:4] <- c("left","right")
queComparisonPol

totalPol <- multi.tab(list(donateBayesLeft,donateBayesRight,infoBayesLeft,infoBayesRight,queBayesLeft,queBayesRight))
names(totalPol)[3:8] <- c("left","right","left","right","left","right")
print(xtable(totalPol),floating=FALSE, include.rownames=FALSE)

#now frequentist
donateFrqLeft <- glmer(modelDonate,
                       data = left,
                       family = "binomial")
donateFrqRight <- glmer(modelDonate,
                        data = right,
                        family = "binomial")
infoFrqLeft <- lm(modelInfo,
                     data = left)
infoFrqRight <- lm(modelInfo,
                      data = right)
queFrqLeft <- lm(modelQue,
                    data = left)
queFrqRight <- lm(modelQue,
                     data = right)
stargazer(donateFrqLeft,donateFrqRight,infoFrqLeft,infoFrqRight,queFrqLeft,queFrqRight,
          star.cutoffs = c(0.05,0.01,0.001), model.numbers = FALSE,
          column.labels =  c("left","right","left","right","left","right"),
          omit.stat = c("adj.rsq", "f","ser","aic","bic","ll"),type="latex")
r2beta(donateFrqLeft)
r2beta(donateFrqRight)
r2beta(infoFrqLeft)
r2beta(infoFrqRight)
r2beta(queFrqLeft)
r2beta(queFrqRight)

##############################
#Frequentist analysis of main text tables
##############################
modelDonate <- donateDum ~ pol * sci + (1|player.orderDonation)

donateFrq <- glmer(modelDonate,
                   data = test,
                   family = "binomial")
summary(donateFrq)
donateFrqWomen <- glmer(modelDonate,
                        data = women,
                        family = "binomial")
donateFrqMen <- glmer(modelDonate,
                      data = men,
                      family = "binomial")
donateFrqYoung <-glmer(modelDonate,
                       data = young,
                       family = "binomial")
donateFrqOld <- glmer(modelDonate,
                      data = old,
                      family = "binomial")
#table 4
library(stargazer)
library(rsq)
library(r2glmm)
stargazer(donateFrq,donateFrqWomen,donateFrqMen,donateFrqYoung,donateFrqOld,
          star.cutoffs = c(0.05,0.01,0.001), model.numbers = FALSE,
          column.labels =  c("Overall","Women","Men","Young","Old"),
          omit.stat = c("adj.rsq", "f","ser","aic","bic","ll"),type="latex")

r2beta(donateFrq)
r2beta(donateFrqWomen)
r2beta(donateFrqMen)
r2beta(donateFrqYoung)
r2beta(donateFrqOld)

modelQuestions <- questions~pol * sci
questionsFrq <- lm(modelQuestions,
                      data = test)
summary(questionsFrq)
questionsFrqWomen <- lm(modelQuestions,
                           data = women)
questionsFrqMen <- lm(modelQuestions,
                         data = men)
questionsFrqYoung <- lm(modelQuestions,
                           data = young)
questionsFrqOld <- lm(modelQuestions,
                         data = old)
#table 5
stargazer(questionsFrq,questionsFrqWomen,questionsFrqMen,questionsFrqYoung,questionsFrqOld,
          star.cutoffs = c(0.05,0.01,0.001), model.numbers = FALSE,
          column.labels =  c("Overall","Women","Men","Young","Old"),
          omit.stat = c("adj.rsq", "f","ser","aic","bic","ll"),type="latex")

r2beta(questionsFrq)
r2beta(questionsFrqWomen)
r2beta(questionsFrqMen)
r2beta(questionsFrqYoung)
r2beta(questionsFrqOld)

questionsComparison

modelInfo <- infoReq ~ pol * sci
infoFrq <- lm(modelInfo,
                 data = test)
infoFrqWomen <- lm(modelInfo,
                      data = women)
infoFrqMen <- lm(modelInfo,
                    data = men)
infoFrqYoung <- lm(modelInfo,
                      data = young)
infoFrqOld <- lm(modelInfo,
                    data = old)
#table 6
stargazer(infoFrq,infoFrqWomen,infoFrqMen,infoFrqYoung,infoFrqOld,
          star.cutoffs = c(0.05,0.01,0.001), model.numbers = FALSE,
          column.labels =  c("Overall","Women","Men","Young","Old"),
          omit.stat = c("adj.rsq", "f","ser","aic","bic","ll"),type="latex")

r2beta(infoFrq)
r2beta(infoFrqWomen)
r2beta(infoFrqMen)
r2beta(infoFrqYoung)
r2beta(infoFrqOld)


